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ABSTRACT 

We develop a detailed chemical model for the starless cores of strongly magnetized 
molecular clouds, with the ambipolar diffusion-driven dynamic evolution of the clouds 
coupled to the chemistry through ion abundances. We concentrate on two representative 
model clouds in this initial study, one with magnetic fields and the other without. The 
model predictions on the peak values and spatial distributions of the column densities 
of CO, CCS, N2H+ and HCO^ are compared with those observationally inferred for the 
well-studied starless core L1544, which is thought to be on the verge of star formation. 
We find that the magnetic model, in which the cloud is magnetically supported for 
several million years before collapsing dynamically, provides a reasonable overall fit to 
the available data on L1544; the fit is significantly worse for the non-magnetic model, in 
which the cloud collapses promptly. The observed large peak column density for N2H"^ 
and clear central depression for CCS favor the magnetically-retarded collapse over the 
free-fall collapse. A relatively high abundance of CCS is found in the magnetic model, 
resulting most likely from an interplay of depletion and late-time hydrocarbon chem- 
istry enhanced by CO depletion. These initial results lend some support to the standard 
picture of dense core formation in strongly magnetized clouds through ambipolar dif- 
fusion. They are at variance with those of Aikawa et al. (2001) who considered a set 
of models somewhat different from ours and preferred one in which the cloud collapses 
more or less freely for LI 544. 

Subject headings: ISM: clouds - ISM: individual (L1544) - ISM: molecules - ISM: mag- 
netic fields - MHD - stars: formation 



1. Introduction 



Through a close interplay between theory and observation, a "standard" picture for the for- 
mation of isolated, low-mass stars has emerged (Shu, Adams &: Lizano 1987). At the heart of this 
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picture lie the so-called "dense cores", studied extensively by Myers and coworkers (e.g., Myers 
1999). These dense cores are intimately associated with star formation, with roughly half of them 
already harboring infrared sources (Bcichman et al. 1986). The other half are termed "starless 
cores" , and they are the focus of our investigation. 

The starless cores are thought to be condensed out of strongly magnetized, turbulent back- 
ground clouds. Direct Zeeman measurements to date, as compiled by Crutcher (1999), suggest 
that the field strength is close to the critical value required for the molecular cloud support, after 
likely geometric corrections (Shu et al. 1999). Core formation can be driven by either turbulence 
decay (Myers 1999) or ambipolar diffusion. We shall focus on the latter, which has been studied 
quantitatively by many authors (e.g., Nakano 1979; Lizano & Shu 1989; Basu & Mouschovias 1994; 
Ciolek & Mouschovias 1994). A general conclusion is that dense cores are formed on a time scale 
several times the free fall time of the background. This relatively long formation time should leave 
a strong imprint on the chemistry of the cores. 

The goal of our investigation is to predict the spatial distributions of various molecular species 
at different stages of cloud evolution. For this initial study, we will concentrate on CO, CS, 
CCS, NH3, N2H+ and 1100+, which are often used to probe the physical conditions of star- 
forming clouds. We take into account the feedback of chemistry on the cloud dynamics, through 
ion abundances, which regulate the ambipolar diffusion and thus cloud contraction rate. Through 
a detailed comparison of the predicted and observed chemical abundances of starless cores, we seek 
to provide further support for the standard picture of isolated star formation involving ambipolar 
diffusion. A starless core well suited for such a purpose is the dense core of L1544, for which an 
extensive data set is becoming available (e.g., Caselli et al. 2001a,b). 

The L1544 starless core is located in the eastern part of the Taurus molecular cloud complex, 
at an estimated distance of 140 pc (Elias 1978). It is an elongated core that shows substantial 
infall motions (up to ~ 0.1 km s~^) on both large (~ 0.1 pc) and small (~ 0.01 pc) scales, based 
on single-disk observations of CO, ^^CO, C^^O, CS, C^^S, HCO+, and N2H+ (Tafalla et al. 1998) 
and interferometric observations of N2H^ (Williams et al. 1999). Intcrferometric observations 
of CCS by Ohashi et al. (1999) reveal a ring-like structure, which shows evidence for collapse 
as well as rotation. The low central emission is most likely due to a depletion of CCS, judging 
from the fact that dust continuum peaks inside the ring (Ward- Thompson, Motte & Andre 1999). 
Indeed, the column density distribution of the L1544 core is rather well determined^, not only 
from the dust continuum emission, but also from absorption against a mid-infrared background 
(Bacmann et al. 2000). The core has a small, high column (and volume) density central plateau 
(of radius ~ 2900 AU), surrounded by an envelope in which the column (and volume) density 
decreases rapidly outward, before it joins rather abruptly with the background at a radius of order 
10^ AU (Bacmann et al. 2000). The plateau-envelope structure has been interpreted in terms 



^See, however, Evans et al. (2001) who find that a singular isothermal sphere cannot be ruled out from the 
modeling of SCUBA data alone, when the lower temperature at the core center relative to edge is taken into account. 
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of the ambipolar diffusion-driven dynamical evolution of a magnetized cloud (Li 1999; Ciolek &: 

Basu 2000); and the large column density contrast between the central plateau and the background 
('^ 20; Bacmann et al. 2000) suggests that the L1544 core is rather evolved (Williams et al. 1999), 
and may be on the verge of forming a star (or stellar system). Further evidence for the core being 
in an advanced stage of evolution comes from IRAM 30m observations of Caselli et al. (1999), who 
found that CO is depleted by a factor of order ~ 10 at the dust continuum peak. These authors also 
noted that the (likely) optically thin lines of D^'^CO"'" and HC^^O"'" are double-peaked, which may 
be another indication that these species are centrally depleted and their emission comes mainly from 
an infalling shell. The wealth of molecular line and dust continuum data, coupled with the recent 
Zeeman measurement of magnetic field strength (Crutcher & Troland 2000) and the determination 
of magnetic field direction from submillimeter polarization measurements (Ward- Thompson et al. 
2000), makes LI 544 a Rosette Stone for theories of dense core condensation leading to isolated 
low-mass star formation. 

Many of the dynamical and chemical features of L1544 are shared by L1498, another quiescent, 
elongated starless core also located in the Taurus molecular complex. Lcmmc et al. (1995) studied 
the core in detail, and found that the distributions of both C^^O and CS are ring-like, indicating 
a central depletion of these two species. The CO depletion factor is inferred to be of order 8 or 
higher (Willacy, Langer &: Velusamy 1998), comparable to that of L1544, based on ISOPHOT 100 
and 200 /xm observations. Single-dish and interferometric observations of CCS (Kuiper, Langer & 
Velusamy 1996) show that its distribution is also ring-like. Moreover, the NH3 emission appears 
to be centrally peaked, again similar to L1544. Perhaps most strikingly, there is evidence for infall 
motions on a large (r^ 0.1 pc) scale in the double-peaked profiles of both optically thick CS lines 
and the likely optically thin C^'^S(2-1) transition (Lemme et al. 1995); the latter could come from 
an infalling shell, as in L1544 (Caselli et al. 1999). The optically thin tracer N2H+ shows, however, 
a single peaked spectrum (Lee, Myers Sz Tafalla 2001), indicating that it is centrally peaked, in 
contrast to C^^S. These similarities to L1544 point to a relatively advanced stage of pre-protostellar 
evolution for L1498, making it another valuable testing ground for chemical and dynamical models 
of (low-mass) star formation. 

Common to both sources is the ring (or shcll)-likc structure of CCS, a well-known "early time" 
species (Bcrgin & Langer 1997) tracing chemically "young" material of order 10^ yrs or less. This 
prominent feature should provide a strong constraint on models. Additional constraints should 
come from the inferred large CO depletion near the core center as well as the spatial extent and 
speed of the infall motions present in both L1544 and L1498. The question we want to address 
is: can the standard picture of dense core formation involving ambipolar diffusion satisfy these 
chemical and kinematic constraints simultaneously? A partial answer to this question has been 
provided by Bergin & Langer (1997), who considered the chemical evolution of a parcel of cloud 
material whose density increases with time in a prescribed way (i.e., a standard time dependent 
but depth independent model). They found that differential depletion of various molecular species 
onto dust grains must be a key ingredient in understanding the differentiated chemical structure 
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observed in L1498. 

Aikawa et al. (2001) went one step further, and followed the evolution of molecular abundances 
of a dynamically collapsing starless core as a whole, with an emphasis on the spatial distribution. 
The study was based on the self-similar solution of Larson (1969) and Penston (1969) and its 
artificially delayed analogs. Among the dynamical models adopted, they found that the undelayed 
Larson-Penston model provides the best overall fit to the observational data on CO, CCS and N2H+ 
for L1544. If this result is robust, its implications would be far reaching. It is well known that 
the Larson-Penston solution describes the collapse of a cloud essentially on a free-fall time scale. 
The fact that it matches observations better than its delayed analogs appears to pose a serious 
challenge to the standard picture of star formation involving strongly magnetized cloud cores, 
which are formed through ambipolar diffusion on a time scale much longer than the free-fall time 
scale. However, the large infall velocity of the Larson-Penston solution, approaching 3.3 times the 
sound speed, is clearly incompatible with those inferred for L1544 and other starless cores, which are 
typically less than the sound speed (Lee et al. 2001). Furthermore, the best-fit model of Aikawa et 
al. under predicts the abundance of N2H''" by a large factor of ~ 20. These discrepancies motivate 
us to reexamine the molecular evolution of starless cores using a more sophisticated dynamical 
model, one that is coupled to the cloud chemistry. 

The dynamical model we will adopt is that of Li (1999). It is a spherical model of starless 
cores that takes into account the dynamic effects of magnetic field approximately. As envisioned 
in standard scenario of low-mass star formation, the cores in the model evolve through ambipolar 
diffusion, and their dynamics are coupled to chemistry through the abundances of charged species. 
Previously, we (Shematovich, Shustov, Sz Wiebe 1997; 1999) have modeled the chemistry of dy- 
namically evolving clouds in detail. The chemical model will be combined with the dynamical 
model into a coupled dynamical-chemical model of magnetized starless cores. An advantage of 
the combined model is that it allows for a determination of not only the time evolution of various 
molecular species but also their spatial distributions and kinematics. We describe the formulation 
of the combined model in §2. In §3, we present representative model results, focusing on the density 
distribution, velocity field and abundance distributions of several commonly used molecular species. 
These results are compared with observations of L1544. We conclude and discuss our main results 
and future refinements in S4. 



2. Formulation of Coupled Dynamical-Chemical Model 

2.1. Dynamical Model 

Ambipolar diffusion-driven evolution of strongly magnetized clouds has been modeled by many 
authors, although a complete understanding is still lacking. The numerical complexity of the 

problem forced most workers to adopt either a quasi-static approximation (e.g., Nakano 1979; 
Lizano &; Shu 1989) or a simplified disk-like geometry (e.g., Basu & Mouschovias 1994; Ciolek & 
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Mouschovias 1994). Some 2D axisymmetric dynamical models are available (Fiedler & Mouschovias 
1993; Dcsch &: Mouschovias 2001), although incorporation of a full-blown chemical network into 
such models would be prohibitively expensive computationally at present. A simpler approach to 
the dynamical problem was pioneered by Safier, McKee & Stahler (1997). By retaining only the 
pressure component of the magnetic forces they were able to retain a spherical geometry, which 
simplifies the dynamics greatly. Indeed, they were able to obtain analytic solutions to both core 
formation and collapse assuming further that the thermal pressure is negligible and that the ions 
(and magnetic field lines tied to them) move much more slowly than neutrals. These last two 
assumptions have been relaxed in Li (1998) and the resulting set of governing equations can be 
solved efficiently, particularly in a Lagrangian form (Li 1999). The formulation of the dynamical 
problem by Li (1999) is employed in our coupled dynamical-chemical model. 

For reference, we list below the equations governing the dynamic evolution of a spherical, 
magnetized cloud (taken from Li 1999) in dimensionless radius (, mass m, density p, velocity u, 
magnetic field strength b and time r: 

^ = -, (1) 

dm pC,"^ ' 

du m a2 ^ P , 



dT e ^ dm V2ac ^ 2 ' ' 
d ( b\ d flAb^C^ db , 

(3) 



\PC J V ^ff P^/^ dm 

The density p and field strength B are scaled by their initial values at the cloud center pc and Be- 
The time t is measured in units of the initial free-fall time 

^^'-^ = (47rGpe)V2' 
at the cloud center, and the velocity V in units of the Alfven speed 

Scales for the radius r and mass M are then Tc = VA,cts,c and Mc = Airpcr^. There are two 
dimensionless quantities in the governing equations. The first, ac, is the ratio of magnetic pressure 
to thermal pressure at the cloud center at t = 0. For simplicity, it will be taken to be unity, 
consistent with the current data on the magnetic fields in the interstellar clouds (Crutcher 1999). 
The second is the magnetic coupling parameter h'Q defined as the ratio of the local free-fall time 
(at a given density p) to the magnetic field-neutral coupling time scale 

^^-kGp 
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where pi is the total mass density of ions, and 7 is an ion-neutral drag coefficient. An approximate 
value for the drag coefficient of 7 = 3.5 x 10^^ cm^ g~^ s~^ is given in Shu (1992), obtained 
assuming a mean molecular weight of ions of ~ 30. We will adopt the above simple prescription 
for the coupling parameter in this initial study, although a more elaborate treatment taking into 
account of contributions from charged dust grains is possible (e.g., Nishi, Nakano k. Umebayashi 
1991; Li 1999) and will be implemented in the future. 



2.2. Chemical Model 

The chemical model used here has been described in detail by Shematovich et al. (1997, 1999). 
Here we summarize the main features briefly. Wc solve the equations of chemical kinetics describing 
reactions in interstellar clouds of relatively low density (~ 10^ — 10^ cm^^) and temperature (T ~ 
10— 100 K). Molecules are formed and destroyed both in gas-phase and on grain surfaces. The dust 
grains are supposed to be well mixed with gas and contain 1% of the cloud mass. We label different 
species, both neutral and charged, by i = 0, 1, . . . , M, with i = denoting electrons. The evolving 
chemical composition is described through the number densities n?(r, t) in gas-phase and (r, t) on 
icy mantles. The exchange between gas and dust is driven by accretion and dcsorption processes. 
In general, the chemical abundances of all species are computed by integrating the equations 

^nf (r, t)=Y,Y. ^§"f - E ^^-"f - ^rnf + fcf ^nf , (8) 

3 I 3 

(r, t) = E E - E + ^r^f - ^""'nt, (9) 

3 I 3 

where Kf- and Kf- are the gas-phase and grain surface chemical reaction rates. The coefficients 
^cc ^-j^^ j^des (jggng i]^Q accretion and desorption rates for the ith component. These equations 
arc integrated, along with the dynamical equations, at each point of the flow. The computed 
abundances are then used to determine the magnetic coupling parameter t'gf , and to obtain column 
densities for various commonly used molecular tracers for comparison with observations. 



2.2.1. Gas-Phase Reactions 

We consider gas-phase reactions in the dense, dark regions of molecular clouds that are well 
shielded from the external UV radiation field. The reactions are initiated by galactic cosmic rays 
which ionize the most abundant species with an adopted rate of 1.3 x 10""^^ s"^. The primary ions 
that are formed react quickly with abundant hydrogen and CO molecules, defining the main paths 
of chemical evolution. 

We use the UMIST Rate95 file (Millar, Farquhar & Willacy 1997). The included elements are 
H, He, C, N, O, S, Si, Na, Mg, and Fe. Our chemical reaction network is mainly limited to those 
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species containing no more than two atoms of trace elements (i.e. all elements other than H and He), 
with a few exceptions. One of these exceptions is C2S, which has been mapped interferometrically 
in L1498 (Kuipcr et al. 1996), L1544 (Ohashi et al. 1999), and others. Its abundance and spatial 
distribution provide a crucial constraint on models, as discussed by Aikawa et al. (2001) and shown 
below. In total, our gas-phase network contains about 2000 reactions between 59 neutral and 83 ion 
species. 

2.2.2. Grain Surface Reactions 

The formation of icy mantles due to accretion and desorption has been taken into account. 
The accretion rate is given by 

= 7ra^ <Vi> SiUgrS'^, (10) 

where a and rigr arc the grain radius (taken to be 10~^ cm) and dust number density, < > is 
the average thermal velocity of the species, and Si the sticking probability. We assume Si = 0.3 
for all neutral species except H2 and He (Willacy, Rawlings & Williams 1994). Two desorption 
mechanisms are included. They are thermal evaporation and cosmic ray induced impulsive thermal 
evaporation (Hasegawa & Herbst 1993). The dust temperature is assumed to be 7 K everywhere, 
close to the value inferred by Evans et al. (2001) and Zucconi, Galli &; Walmsley (2001) near the 
center of L1544. The values of the sticking probability and desorption energies are uncertain. The 
sensitivity of our model to these uncertainties will be explored elsewhere. 

After being accreted onto the grain surface, the species are allowed to react with one an- 
other mainly through the hydrogen addition reactions, with rates taken from Hasegawa, Herbst &; 
Leung (1992). The surface products include saturated molecules (such as CH4, H2O and NH3), 
homonuclear molecules (such as C2, N2, and O2) and their hydrogenated forms, among others. The 
accretion limit to the surface hydrogen addition chemistry was taken into account through the first 
modification of rate equations proposed by Caselli, Hasegawa, & Herbst (1998). In addition to 53 
surface reactions (inchiding molecular hydrogen formation), we also consider the recombination of 
ions on grain surfaces. The products of dissociative recombination are assumed to return into the 
gas phase immediately. The effects of surface reactions on abundance distributions will be fully 
examined in a subsequent paper. 

2.3. Initial and Boundary Conditions 

We adopt as a starting point for our calculations an isolated molecular clump in mechanical 
equilibrium, with a magnetic pressure equal to the thermal pressure everywhere. The gas tempera- 
ture is assumed to be 10 K and is kept constant throughout the computations. The initial molecular 
hydrogen central density is set to nu^ = 10^ cm~^ for all models. Without ambipolar diffusion, the 
cloud would remain in the equilibrium forever. At i = ambipolar diffusion is switched on. The 
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cloud evolution is then followed numerically with reflection boundary conditions (i.e., zero spatial 
gradient) at the cloud center and with a free pressure boundary at the cloud outer edge. 

The initial chemical composition is assumed to be mostly atomic. Species with non-zero initial 
abundances are listed in Table 1. Their abundances are given with respect to the number density 
of hydrogen nuclei. All hydrogen is assumed to be in a molecular form initially. 

2.4. Numerical Method 

We divide the computational domain into spatial zones of equal mass. At each time step, we 
first advance the chemical model to evaluate ion density p\ in each spatial zone. Equations (1), 
(2), and (4) are then solved to update the velocity, radius of zone boundary and density. A new 
distribution of magnetic field is obtained by solving the diflFusion equation (3), with the magnetic 
coupling parameter calculated from the ion densities at the beginning of the time step. 

The value of time step is computed from the usual magnetic Courant-Priedrich-Lewy condition 
to ensure dynamic stability. The chemical model is advanced over this time interval independently, 
with its own algorithm and smaller time substeps that are determined from the required abundance 
accuracy. Wc solve the chemical kinetics equations with the standard LSODE package, and termi- 
nate the computation when the central density reaches 71^2 = 10^ cm~^. This density is close to 
the value inferred at the center of LI 544. At higher densities the effects of charged dust grains on 
may become dominant (Nishi et al. 1991), and such effects have not been taken into account in 
the current version of the code. 

3. Representative Model Results 

The initial conditions for forming starless cores are not well constrained, either observationally 
or theoretically. There is in principle a large parameter space to be explored. In this initial study, we 
will focus on two models that highlight the effects of magnetic field on the dynamics and chemistry. 
In the first model, we consider a standard, magnetic cloud with a mass of 20 Mq (StM+B). The 
second model is identical to the first except that the magnetic field strength is set to zero (i.e., 
non- magnetic; StM-B). For reference, we also consider the dynamics of a model that is the same 
as the first except that the magnetic coupling parameter is fixed at the canonical value of = 10 
instead of being computed self-consistently. These models are summarized in Table 2. 

3.1. Dynamics 

The main dynamic quantities of the models (density distribution and velocity field) are dis- 
played in Figs. 1-3, along with the ionization fraction and magnetic coupling parameter for the 
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StM+B model. As mentioned earlier, all models start with a central density of tzhz = 10^ cm~^, 
typical of clouds in transition from being diffuse to dark. Note that when a certain center den- 
sity is reached, the difference in density profile is modest, with that in the non-magnetic model 
somewhat peakier than those in the magnetic models. The difference in the time it takes to reach 
a given density is more pronounced. In Table 3 we list the time intervals taken to increase the 
central density by successive factors of 10. The time for the first factor-of-10 increase, to a value 
of 10*^ cm~^ characteristic of dense molecular cores, is ~ 1 Myr for the non-magnetic model. It is 
substantially shorter than those for the magnetic models ('^ 4 — 5 Myrs). The times for subsequent 
10-fold increases become shorter and shorter, and the difi"crcnces in time between the models also 
get smaller. The time for the second (and third) factor-of-10 increase is longer for the magnetic 
models than for the non-magnetic one by only a factor of ~ 2.5 (and ~ 2). In other words, after 
formation a magnetic starless dense core typically lasts for only a few free-fall times before stellar 
birth. That the evolutionary time scale becomes closer to that of a non-magnetic, free-falling model 
at higher densities is a,n im,portant feature of the ambipolar diffusion- driven cloud evolution. This 
characteristic cannot be captured by simply adding a constant delay factor to a non-magnetic m.odcl, 
as done for example in Aikawa et al. (2001). The differences in time scale among the models will 
have a decisive influence on the chemistry, with potentially observable consequences. 

Another potentially observable difference among the models is the velocity field. From Figs. 1- 
3 we find that the infall speed is typically a few tenths of 1 km s~^, comparable to the isothermal 
sound speed of the cloud (a = 0.188 km s~^ at 10 K). As expected, the non-magnetic model 
collapses right away, approaching a maximum speed ~ 2.3 times larger than the sound speed a by 
the time of a lO'^-fold increase in the central density. Even this maximum speed is substantially 
lower than that in a Larson-Penston flow, which is 3.3 times the sound speed. The infall speed 
in magnetized models increases more gradually, reaching a maximum value of ~ 0.35 km s~^ (or 
about twice the sound speed) for the StM-|-B and reference model when the central density reaches 

= 10^ cm~^, roughly the density inferred for the central region of L1544 (Ward-Thompson et 
al. 1999). Note that the density distribution flattens near the center for all three models, more so 
in the magnetic models than in the non-magnetic one. The radius of the central flat region is of 
order ~ 2 — 3 x 10^ AU for the StM-|-B model, consistent with those deduced by Ward-Thompson 
et al. (1999) and Bacmann et al. (2000). The predicted peak H2 column density of ~ 10^^ cm"^ 
also agrees with the value inferred by Ward-Thompson et al. (1999) within a radius of ~ 900 AU. 
However, the observationally inferred infall velocity for L1544 of order 0.1 km s~^ (Tafalla et al. 
1998; Williams et al. 1999; Ohashi et al. 1999) is lower than predicted. The discrepancy, although 
substantial, is less severe in our model than in those of Aikawa et al. (2001). It could be further 
reduced by taking into account the projection effects in a more realistic flattened geometry and the 
effects of magnetic tension on cloud dynamics (Ciolek k, Basu 2000). 
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Table 1: Species with Non-Zero Initial Abundances 



Species Abundance 



H2 


0.5 




He 


0.1 




Na 


2.1( 


-6) 


Mg 


1.5( 


-7) 


Fe 


1.0( 


-7) 


C 


1.0( 


-5) 


N 


2.2( 


-5) 





1.3( 


-4) 


CO 


4.0( 


-5) 


s 


8.0( 


-8) 


Si 


8.0( 


-9) 



Tal")k^ 2: Model SumuiMry 


Model Name 


us 


Magnetic Field 


Reference 


fixed at 10 


yes 


StM+B 


variable 


yes 


StM-B 


not relevant 


no 



Table 3: Time Intervals (in Myrs) for Central Density Increase 



H2 density (cm ^) 


Reference Model 


StM+B 


StM-B 


10^-10^ 


4.35 


5.05 


1.32 


lO'^-lO^ 


0.47 


0.53 


0.19 


10^-10^ 


0.10 


0.11 


0.05 
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3.2. Ionization Fraction and Magnetic Coupling Parameter 

Before showing results on chemical abundances, we comment on the ionization fraction Xi 
(defined as the ratio of the number density of all ions to that of hydrogen nuclei) and magnetic 
coupling parameter ug, both of which are shown in Fig. 2 for the StM+B model. The fraction Xi 
decreases with density, as expected. Its lowest value of ~ 3 x 10^^ at the center of the cloud is 
somewhat higher than that inferred for L1544 by Caselli et al. (2001b), although the inference is 
model dependent. The ionization fraction yields a magnetic coupling parameter vg within a factor 
of two of the canonical value of 10. 



3.3. Chemical Abundances and Comparison with L1544 

3.3.1. Observationally Inferred Abundances in L1544 

We concentrate on the abundances of six species, including CO, CS, CCS, HCO"*", NH3 and 
N2H+. Other species will be considered in a subsequent paper of the scries. For CO, CCS, HCO^, 
and N2H+, the column density distributions have been inferred from observations in L1544. Caselli 
et al. (2001b) found that CO is heavily depleted near the center of L1544, by a factor up to ~ 10. 
Its column density distribution is flat or slightly depressed near the center, with an average central 
value of ~ 1.5 x 10^^ cm~^, although large scatter exists in the data (see their Fig. 5). Both 
HCO"*" and N2H"'' arc centrally peaked, with a peak column density of ~ 10^^ and ~ 2 x 10^^ cm~^ 
respectively. Ohashi et al. (1999) mapped L1544 intcrfcrometrically in CCS. They found that CCS 
peaks in a ring of radius ~ 7, 500 AU. The peak CCS column density is ~ 4 x 10^^ cm~^, which 
is ~ 1.4 times higher than the value at the center (Aikawa et al. 2001). Observational data on 
NH3 and CS appear to be less complete. Using the Nobeyama 45 m telescope, Suzuki et al. (1992) 
observed L1544 in NH3 and derived an average column density of 2.7 x lO^^cm""^ within a region of 
radius ~ 6, 000 AU, corresponding to a beamwidth of 80". L1544 has been observed in CS (Tafalla 
et al. 1998), although the peak value and spatial distribution of its column density are not yet 
available. 



3.3.2. Standard Magnetic Model 

The column density distributions of the six species of interest for the StM+B model, obtained 
by integrating the number densities along the line of sight at various distances from line of sight 
through the cloud center, are shown in Fig. 4. We will concentrate on the last time (t=5.69 Myrs) 
when the density distribution fits approximately that observed (Ward-Thompson et al. 1999). It 

is immediately evident from the figure that a central hole exists in the distribution of CCS, with 
a radius ('^ 7,000 AU) and depth (~ 1.3) that match almost exactly those observed (Ohashi et 
al. 1999). Note the differentiated chemical structure, with the centrally peaked NH3 and N2H"'" 
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surrounded by a CCS ring, which is in turn embedded within a sHghtly broader CS ring^. The 
predicted peak cohimn density for N2H"'' of 1.3 x 10^'^ cm^^ is within a factor of 2 of the observed 
value (~ 2 X 10^"^ cm^^). The predicted peak column density for NH3 of 5.5 x 10"*^^ cm^^ is about 
twice the value 2.7 x 10^"^ cm~^ given by Suzuki et al. (1992). Averaging the column density over a 
radius of 6,000 AU (the beamwidth) yields a value ~ 2 x 10^^ cm~^, which is in a closer agreement 
with Suzuki et al.'s value. Prom Fig. 4 we see that HCO"*" is also centrally peaked, although less 
sharply than N2H"'", again in agreement with observation (see Fig. 5 of Caselli et al. 2001b). Its 
predicted peak column density of 4.5 x 10^^ cm~^ is about a factor of two lower than observed 
(~ 10^^ cm^^). The CO column density in this model shows a relatively flat distribution near the 
center, with a maximum value of 1.6 x 10^^ cm~^, almost identical to the observed value. Its spatial 
distribution appears to be consistent with that inferred by Caselli et al. (2001b) within errors. The 
model predicts that the column density distribution of CS should have a shallow depression near 
the center. Data for testing this prediction are not yet available. We conclude that, within a factor 
of two or so, the predictions of the StM+B model on the abundances of CO, N2H+, HCO+ and 
NH3 are consistent with the data available on L1544. 

The abundance and spatial distribution of CCS deserve special attention. We have mentioned 
the close fit to the observed spatial distribution and will now concentrate on abundance. The CCS 
molecule is a well-known "early time" species. Its abundance should decrease, according to Bergin 
& Langer (1997), rapidly after a few times 10^ years, as a result of neutral atomic carbon being 
converted into CO (see also Suzuki et al. 1992 and Aikawa et al. 2001). One might expect to 
find little CCS in the magnetic cloud of the StM+B model, which reaches the observed state after 
nearly 6 million years of evolution. However, this turns out not to be the case, as shown in Fig. 4. 
The predicted peak CCS column density is 1.3 x 10^^ cm^^, about a third of the observed value 
(~ 4 X 10^^ cm~^). The agreement would be closer, if one takes into account of likely projection 
effects. The CCS ring is seen almost edge-on, judging from its large aspect ratio (~ 3; Ohashi et 
al. 1999). The nearly edge-on geometry implies a smaller peak column density perpendicular to 
the ring, by a factor of ~ 3. It turns out that CCS plays a crucial role in discriminating models. 
We will postpone a more detailed discussion of this molecule to § 4.1. Here, we simply note that 
the StM+B model appears to reproduce all abundance data on L1544 reasonably well. It is the 
"standard" model against which other models will be compared. 



3.3.3. Free-Falling Non-Magnetic Model 

Without magnetic support, the cloud considered in the previous subsection collapses promptly, 
reaching the observed state in 1.56 million years. This faster dynamic evolution affects the chemical 



^This is not exactly what is observed in L1498, where CS appears to concentrate near the inner edge of the CCS 
ring (Kuiper et al. 1996). The situation with L1544 is unclear, although there is some indication that a central hole 
exists in CS as well (Tafalla 2001; priv. comm.). 
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abundances, as shown in Fig. 5. The effects on CO, NH3, N2H"'" and HCO"*" are relatively modest. 
The peak column density of CO differs little from that of the standard model. Its spatial distribution 
appears more centrally peaked, however, which is more difficult to reconcile with observation. The 
distributions of N2H+, NH3 and HCO"*" are centrally peaked, as in the standard model. Compared 
with the standard model, the peak value for N2H+ is down by a factor of ~ 2, making it significantly 
below the observed value (by a factor of ~ 4) . The column density of NH3 is also down by a similar 
factor of ~ 2, making the fit to the value inferred by Suzuki et al. (1992) worse. The column 
density of HCO"*" is lower by a modest factor of ^ 20%. 

The difference in CS and CCS is more pronounced. Unlike in the standard model, the column 
density distributions of CS and CCS in the non-magnetic model arc centrally peaked instead of 
showing a central depression. While the observational situation with CS is still unclear at present, a 
central peak in CCS is definitively ruled out. One potentially desirable feature of the non-magnetic 
model is that the peak column density of CCS is about twice the value in the standard model, 
making it closer to the observationally inferred value without projection corrections. However, 
its overall fit to all available data is significantly worse than that of the standard model. This 
result lends some support the standard picture of low-mass star formation involving core formation 
through ambipolar diffusion over a time scale substantially longer than free fall. The peak value 
and spatial distribution of the column density of CS would further constrain the models but are 
not yet available. They are urgently needed. 

We note that our non-magnetic model is similar to the undelayed model of Aikawa et al. (2001), 
except that we have included some surface reactions and adopted a smaller sticking probability. 

4. Discussion and Conclusions 

4.1. Magnetic Cloud Support and CCS Abundance 

The success of the standard, magnetic model in fitting the data on L1544 hinges to a large 
extent on the relatively high abundance predicted for CCS. This result is somewhat surprising in 
view of the fact that CCS is a well-known "early time" species, whose abundance should decline 
rapidly after a few times 10^ years according to Suzuki et al. (1992), Bergin & Langer (1997), and 
Aikawa et al. 2001). Part of the reason for the persistently high CCS abundance in the standard 
model even after some 5.69 Myrs of evolution (the time it takes the cloud to reach the observed 
state) is that we adopted an initial cloud density of nn^ = 10^ cm^"^, which is lower than those of 
Bergin & Langer (10^'^ cm~^) or Suzuki et al. and Aikawa et al. (10"^ cm~^). As a result, the bulk 
of the 5.69 Myrs is spent in the "pre-dense core" phase when the central density remains below 
10^ cm~^ (see Table 3). Only a small fraction (~ 10%) of that time is spent at densities above 
10^ cm~^. That a cloud spends a relatively short time (a few times the free-fall time) at high 
densities after forming a magnetically supercritical core is a general feature of ambipolar diffusion- 
driven evolution. It helps to keep the CCS abundance higher than one would expect based on the 
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total time of cloud evolution. 

The fractional abundances of CCS are shown in Fig. 6 as a function of radius for both the 
magnetic and non-magnetic models. Note the "humps" on the distributions of CCS abundance in 
the magnetic model, which are not apparent in the non-magnetic model. The exact origin of the 
humps is unclear. We suspect that they are related to the depiction of CO in the central high 
density region, similar to the "depletion" peak proposed by Ruffle et al. (1997) for HC3N. Ruffle et 
al. (1997, 1999) showed that C2H and HC3N in dense, cold cores are characterized by a secondary 
"late-time" maximum in their fractional abundances. The second peak is thought to be mainly 
caused by the increased depletion of CO, which allows 0+ ion to react more readily with H2 than 
with oxygen bearing species. The net result is an increase in the rates of production of CH and 
other carbon-bearing molecules (without oxygen). It was found that high late-time C2H and HC3N 
maxima are achieved only when the freeze-out time scale is long compared to the chemical time 
scale. Because C2H and other (neutral and ionized) late-time hydrocarbons are the precursors of 
CCS formation in the reactions with neutral and ionized sulfur through 

C2H, C2H2, C2H3, ... + S+ ^ HCCS+ + (11) 

C2H+ C2H+ C2H+ ... + S ^ HCCS+ + (12) 

HCCS+ + e ^ CCS + H, (13) 

CCS should behave in a way similar to the C2H and HC3N. The competition between depletion and 
late-time hydrocarbon chemistry enhanced by CO-depletion appears to be the most likely cause of 
the humps on the radial profiles of CCS abundance. 

In any case, the humps have apparently kept the CCS abundance of the magnetic model close 
to, and in some regions exceeding, that of the non-magnetic model. They are largely responsible 

for the higher-than-cxpected value of the CCS column density in the magnetic model. In addition, 
the decline of CCS abundance outside the hump, coupled with a steep decrease in the hydrogen 
number density with radius, may explain the observed rapid decline of the CCS column density 
outside the ring (Ohashi 2001; priv. comm.). 

4.2. Depletion of Molecules onto Dust Grains 

4.2.1. Depletion of CO and Other Species 

In the standard model, CO is heavily depleted in the central high density region. The depletion 
is most clearly seen in the first panel of Fig. 7, where the fractional abundance is plotted as a function 
of radius. The abundance is significantly below the canonical value of 4.7 x 10~^ inside a radius of 
~ 10^ AU, by a factor up to ~ 30. The size of the depletion region is roughly consistent with that 
deduced by Caselli et al. (2001b; their Fig. 2). The depletion is less evident in the CO column 
density distribution (see Fig. 4), which is fiat near the center. The ratio of CO and hydrogen in 
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peak column density is a factor of ~ 5 below the canonical value. This depletion factor is about half 
of the value ('^ 9) inferred by Caselli ct al. (2001b). Their inference is based on the observations 
of C^^O, C^^O, and dust emission (Ward-Thompson et al. 1999). In particular, an dust opacity of 
0.005 cm^g~^ was adopted to convert the dust continuum flux into a hydrogen column density. At 
densities as high as 10^ cm~^, an opacity of ~ 0.01 cm^g"-*^ may be more appropriate (Ossenkopf 
&; Henning 1994), which would lower the estimate of the hydrogen column density by half and 
bring a closer agreement between the model prediction and observation. Alternatively, a flattened 
geometry, as observed for L1544, and/or a somewhat higher CO adsorption energy (see Aikawa et 
al. 2001) may enhance the CO depletion in the column density distribution. 

Besides CO, there are other species that are strongly depleted in the high density central 
region. These include CS, CCS and to a lesser extent HCO"*", as can be seen from Fig. 7. The 
nitrogen bearing species, N2H"'" and NH3, are on the other hand hardly depleted; if anything, 
their abundances are slightly enhanced near the center. This differential depletion pattern is in 
agreement with previous findings (e.g., Bergin &; Langer 1997; Aikawa et al. 2001). It may have 
profound effects on the molecular line profiles used to probe cloud kinematics. 

4-2.2. Effects of Depletion on Line Profiles 

An interesting observational fact about L1544 is that some optically thin lines in L1544, such 
as CCS (Ohashi et al. 1999) and HCi^O+ (CaselU et al. 2001b), are double-peaked. As noted by 
Caselli et al. (1999, 2001b), the double-peaked optically thin lines can be explained in a collapsing 
cloud provided that the line-emitting molecules are strongly depleted near the center. We have 
seen from Fig. 7 that CCS and, to a lesser extent, HCO"*" (and thus HC^^O"^) are indeed centrally 
depleted in our standard model. A unique strength of the coupled dynamical-chemical model is 
that it allows for a simultaneous determination of the spatial distribution of a given species and the 
velocity field, the two ingredients for line profile modeling. In a subsequent paper, we will produce 
synthetic line profiles for both optically thin and thick lines (such as CS and HCO^; Tafalla et al. 
1998), which should enable us to constrain the models of core formation leading to star formation 
using both the spatial and kinematic information of molecular species. 

4.3. Conclusions 

We have developed a coupled dynamical and chemical model for the starless cores of strongly 
magnetized molecular clouds. The coupling is achieved through ions, the abundances of which 

are determined sclf-consistently from a chemical network. The ionic abundances control the time 
scale of cloud dynamic evolution via ambipolar diffusion. We have concentrated on a representative 
model in which an isolated magnetic cloud increases its central number density by a factor of 10^ 
to = 10^ cm^^ in 5.69 million years, with an eye on explaining the observational data on the 
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chemical abundances of the well-studied starless core L1544. We find that the predicted peak values 
of CO, CCS, N2H"'" and HCO"*" column densities are within a factor of two or so of those inferred for 
L1544 from observations. The spatial distributions of these species are also consistent with those 
observed. An alternative, non-magnetic model was also considered for comparison. With the same 
initial conditions as in the magnetic model, the non-magnetic cloud collapses promptly, reaching 
the observed state in merely 1.56 million years. It produces a worse overall fit to the available data 
on L1544. In particular, the column density distribution of CCS predicted by the model is centrally 
peaked, which is not observed. 

We conclude that our initial results of modeling lend some support to the standard picture 
of dense core formation out of strongly magnetized molecular clouds involving ambipolar diffusion 
over several dynamic times. There are, however, a number of uncertainties to which the results 
may be sensitive to, including the sticking probability and adsorption energies for various species. 
A parameter survey is needed to firm up the conclusion. 

4.4. Future Refinements 

There are several aspects of the coupled model that we wish to improve upon in the future. 
These include both dynamics and chemistry. Even though our spherical model captures the essence 
of the ambipolar diffusion-driven cloud evolution, the effects of magnetic tension (which tend to 

flatten a cloud) cannot be treated. Indeed, the opposite extreme, a disk-like geometry, may describe 
the observed elongated mass distribution in L1544 better. It should be straightforward to extend 
the model to this geometry, which would allow for the inclusion of rotation whose presence has 
been inferred for L1544 by Ohashi et al. (1999) based on CCS data. On the chemistry side, a 
major uncertainty is the adsorption energies for CO and other species. Different sets of adsorption 
energies have been considered by Aikawa et al. (2001), and a similar parameter study is needed 
for our model. Also uncertain are the probability of molecules sticking onto dust grains and grain 
surface processes, and an exploration of different cases is desirable. Another area of improvement 
would be the treatment of magnetic coupling coefficient, by including the effects of charged dust 
grains (especially small ones) and external UV radiation field. 

The authors are grateful to P. Murphy and B. Turner for access to NRAO computing facilities, 
and to P. Caselli and the referee for helpful comments. This research is supported in part by RFBR 
Grant 01-02-16206. 
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Fig. 1. — Radial profiles of the number density of hydrogen nuclei (a) and infall velocity (b) for the 
reference model at three times when the central density reaches 10, 10^ and 10^ times its initial 
value respectively. 
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Fig. 2. — Radial profiles of the number density of hydrogen nuclei (a), infall velocity (b), ionization 
fraction (c), and magnetic coupling parameter (d) for the standard, magnetic model. 
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Fig. 3. — Radial profiles of the number density of hydrogen nuclei (a) and infall velocity (b) for the 
non-magnetic model. 
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Fig. 4.— Column density distributions of CO, CCS, CS, NH3, N2H+, and HCO+ for the standard, 
magnetic model. 
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Fig. 5.— Column density distributions of CO, CCS, CS, NH3, N2H+, and HCO+ for the non- 
magnetic model. 
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Fig. 6. — Radial profiles of the fractional abundances of CCS for both the magnetic (heavy lines) 
and non-magnetic (light lines) models. 
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Fig. 7.— Radial profiles of the fractional abundances of CO, CCS, CS, NH3, N2H+, and HCO+ 
for the standard, magnetic model. 



